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Abstract. Explicit solution for the 2-point correlation function in a non-equilibrium 
steady state of a nearly isotropic boundary-driven open XY spin 1/2 chain in the 
Lindblad formulation is provided. A non-equilibrium quantum phase transition from 
jrt . exponentially decaying correlations to long-range order is discussed analytically. In the 

tr. ' regime of long-range order a new phenomenon of correlation resonances is reported, 

where the correlation response of the system is unusually high for certain discrete 
values of the external bulk parameter, e.g. the magnetic field. 
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1. Introduction 

In recent years we are witnessing an increasing activity in non-equilibrium physics of 
interacting many-body problems. The reason may be at least two fold. Namely, on 
one hand this type of problems has nowadays become amenable to a detailed real 
lab experiments, e.g. in the context of cold atoms and optical lattices [1], and on 
the other hand, they are arguably connected to important open challenges in solid 
state physics, such as solid-state quantum computation or high temperature super- 
conductivity. Essentially, there are two basic approaches to mathematical treatment of 
non-equilibrium many-body-physics. Perhaps more common approach [2] is to consider 
a large (ideally infinite) system in an equilibrium state, then at some instant of time 
perform a quench of the Hamiltonian or join (couple) several infinite pieces of the system 
in distinct equilibria, and observe what happens in the course of time. For example, 
there might exist stationary long time behaviour. The other, more explicit approach 
[31 H], which we consider here, is to couple a finite (though perhaps large enough) 
system to several external reservoirs which we may describe effectively in terms of a 
master equation (by tracing out their explicit degrees of freedom), and consider a non- 
equilibrium steady state (NESS) to which a central system converges after a long time. 

There is a variety of techniques involving several levels of assumptions and 
approximations in deriving the effective system's dynamics after tracing out the 
reservoirs degrees of freedom [3l H], which at the end result in a simple local- in-time 
(Markovian) linear differential equation for the system's density matrix, the so-called 
quantum Liouville equation. The more general Markovian form of such equations 
is sometimes referred to as the Redfield equation, whereas the more mathematically 
appealing form which manifestly conserves the positivity of the density matrix (and 
can be derived from the Redfield model with an additional, the so called secular or 
rotating wave approximaiton, is the famous Lindblad equation [5]. It has been noted that 
Lindblad driven quantum chains provide a fruitful ground for studying non-equilibrium 
phase transitions [6]. 

Despite the fact that many simple cases of the Redfield and Lindblad master 
equations have been studied, a very few exactly solvable instances where the central 
system consists of many interacting particles have been discovered so far. In particular, 
in the context of addressing the quantum transport problem [Tj |8], one is interested 
in the quantum chain which is coupled to effective (thermal, chemical, or magnetic) 
reservoirs only at the ends of the chain (i.e. via the first and the last spin/particle). 

As an important simple but notable class of exactly solvable quantum Liouville 
equations one can identify open XY spin chains, or open quasi-free fermionic systems 
(to which XY chains can be mapped via Jordan- Wigner transformation) in general. As 
it has been shown in p], the quantum Liouvillean of a general quadratic system of n 
fermions can be explicitly diagonalized using the generalized (non-unitary) Bogoliubov- 
like transformation, and all the properties of NESS and the relaxation process can 
be computed explicitly in terms of diagonalization of a 4n x 4n matrix. The main 
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conceptual tool in this approach is the introduction of the Fock-space structure over 
the space of operators which the density matrix is a member of. This method has been 
later generalized also to more general Redfield master equations [TO] and put to a more 
rigorous mathematical footing [11]. It is remarkable that explicit solution for the NESS 
density matrix can be found also in some cases which go beyond quasi-free models, e.g. 
in XX spin 1/2 chain with on-site dephasing noise [12] or in the XXZ spin 1/2 chain 
for the extremal (strong) driving [13] . On the side of numerical methods, the Liouville 
space version of the time-dependent density matrix renormalization group [TU [I5] was 
successfully implemented to simulate steady states of open quantum chains [16]. For 
some other related explicit results on open quantum chains see e.g. Refs. [T7] . 

As it has been shown in [18], an open XY chain (either in the local Lindblad or 
Redfield pilj setting) exhibits a quantum phase transition from exponentially decaying 
correlations to long range order in NESS when the (bulk) parameters of the model 
are varied. Although a simple heuristic picture in terms of dispersion relation of the 
quasi-particle normal modes has been proposed to explain the transition, the precise 
mathematical and physical understanding remained lacking. In this paper we show that 
the non-equilibrium transition can be explained with explicit analytical calculation in 
the special regime of small anisotropy, where the perturbation theory in the anisotropy 
parameter can be successfully apphed. We note that the small anisotropy open XY 
chain has been recently also studied numerically in terms of the Keldysh formalism [19] , 
where the non-equilibrium phase transition has been re-confirmed. 

Even though being surprisingly simple, our results analytically reproduce all the 
features of the transition to long range order under the additional assumption of 
weak coupling to the baths: namely we (i) give explicit expression for the 2-point 
correlation function, (ii) evaluate the decay rate of the correlation function in the short- 
range regime, and (iii) predict a new phenomenon in the long-range regime, the so 
called resonances of the correlation function at particular discrete values of the bulk 
parameter, e.g. the transverse magnetic field. At these resonant values of the magnetic 
field, in the long range order regime, the response of the system to external driving is 
particularly large as characterized by spontaneous emergence of strong correlations over 
large distances. 

The paper is organized as follows. In section 2 we review the main concepts of 
quantization in the Fock space of operators, write down explicit equations of motion, 
and derive the so-called Lyapunov equation, as the dynamical equation for the 2-point 
correlation function which completely determines NESS. We stress that the results of this 
section represent certain simplification and improvements over the formulation of the 
previous papers [9l [10] , namely we write for the first time a uniform quadratic form for 
the Liouvillean which remains valid in both, even and odd parity sectors of the operator 
space. This we do with a small trick, a simple redefinition of the canonical Majorana 
maps over the operator space. In section 3 we then solve perturbatively, uniformly 
in two small parameters, the anisotropy and the system-bath coupling strength, the 
Lyapunov equation for the correlation (or covariance) matrix for the nearly isotropic 
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open XY chain. We finally discuss our results and conclude in section 4. 

2. Equation of motion in operator Fock space 

In this section we make a very brief review of the technique of Liouvillean diagonalization 
of quadraric fermi systems in the Fock space of operators |9l HD] with an emphasis 
on the equation of motion and real-time djTiamics. We treat a finite system with n 
fermionic degrees of freedom, described by 2n anti-commuting hermitian operators Wj 
(j = 1, 2, . . . 2n). The Redfield master equation (see e.g. [4j) Ij 

^ = £p:=-(iad/7)p + Pp (1) 

serves us as the most general (Markovian) equation of motion that we are able to 
treat [TU] (quantum Liouville equation). The first term is anti-hermitian — (i adH)p := 
—i[H, p] and generates the usual unitary von Neumann equation. The second term - the 
dissipator map - has a memoryless kernel with the following general form [3] 

/■oo 

v = Y, dtr^^^(t)[x^(-t)p,x,] +h.c., (2) 

where F^^ (t) denotes the environment (bath) correlation function and X^{t) denotes 

the Heisenberg picture of hermitian coupling operators X^ (following notation of [TO]). 

The Hamiltonian can be a general quadratic form in Wj and the coupling operators 
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The 2n x 2n matrix H is antisymmetric. Throughout this paper x = (a;i,X2, . . .)""" will 
designate a vector (column) of appropriate scalar valued or operator valued symbols x^- 

2.1. Bilinear form of the complete Liouvillean 

Diagonalization of the Liouvillean C is done in the 4" dimensional Liouville space of 
operators /C, which has the structure of the Fock space with the canonical basis 

p^ = 2-^i^wTwT . . . wZ\ oij e {0, 1}. (4) 

We will use the bra-ket notation for the operator space /C with the inner product 

{x\y) = ti{x^y). (5) 

Further we define the creation/annihilation maps by cl\Pa) = (1 — aj)\wjPa), Ci\Pa) = 
aj\wjPa), which satisfy canonical anti- commutation relations (CAR), {ci,cA = Sij, 

I We use units in which Planck's constant h= 1. 
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{c|,c^} = {cj,Cj} = 0. The Redfield dissipator can be written as a quadratic form in 
these maps 

2n 

^ = E E ^'^.'^ [''^A^k + 4A,k) , (6) 

V k,j=l 

where z^j are complex constants defined as 

/■oo 

z,:=J2 dtexp(-4iHt)x/J^(t) (7) 

and 

^'j,k\x) ■=\[wjX,Wk]), C'lklx) ■.= \[wk,xwj]) (8) 

are fundamental basis dissipators which evaluate to 

£' =(i + v) ic]cl - clc,) +(l-v) {cfcu - cue]), (9) 



C'l, = (i + P) idle] - die,) + (i - P) (CkC, - Ckc]). 

Note that the parity map V = exp (ivr X]j=i c]cj j and the positive and negative parity 

projectors P^ = ^(I ± V) commute with the Liouvillean, [£, P*^^)] = 0. Plugging 
the definitions @ in the equation for the Redfield map (|6]) we see that the dissipator 
decomposes into a positive {'D^) and a negative {'D^) parity components 

V = -1)+^+ + v-p- ^ (10) 
V+ = 2^ ■ (M^ + M*) c^ - 2c^ • (M + M*) c, 

i)- = 2c- (M^ + M*) c - 2c ■ (M + M*) c^ 

where M is a hath-matrix which, can be compactly written as M := Yliv^v®^u^ ^^^ M* 
and M""" denote respectively, the complex conjugate and transpose of the matrix. The 
Liouvillean can be written in a convenient form using An fermionic Majorana maps 

«ij := -7|(cj + c]), alj := — =(cj - c]) and a = {al,al) (11) 

satisfying CAR {a°^j,a,lk} = S^i,uSj,k- Straightforward calculation shows that the 
dissipator has the following form 

V^ = a^- Ai lOi + 0,2 ■ Ag^sfe + Oi ■ A^2Ql2 + Ql2 " ^2,iQli - ^ol, (12) 

V = Oi ■ Ai iCti + 0,2 ■ P^2,2Q^2 - QlI ■ Ai^sfe - fe ■ ^2,\Ql\ " ^0^, 

where we used the definitions 

A°i := M^ - M, A°2 :=iM'^ + iM*, (13) 

A°i := -iM'f - iM, A°2 := M^ - M*, 
Ao:=tr(M + M*). 

The unitary part of the Liouvillean is also quadratic in the Majorana maps a^, Sg [H [10] 

- i adff = -4ic^ ■ He = -2i(a° • Ha° + a° ■ Ha"). (14) 
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Collecting the results (ITU]) . ( 1T^ the complete Liouvillean decomposes in the same way 
as its dissipative part 

C =£+P+ + £-p-, (15) 

= ai- Ai^io-i + 0-2 ■ A2,2ai2 + % ■ Ai^2(h. + fe ' ^-2,10:1 - AqI, 
C^ = a° ■ Ai^io,^ + 0-2 • A2,2ai2 " Qli ■ ^1^202 — Og • A2,iai — AqI, 
where we introduce 2n x 2n matrices 

A^,, := -2i5^,,H + A°,, fx,u = l,2. (16) 

Summing the odd and even parity Liouvilleans (TT5l) and using the definition of the parity 
projectors we obtain 

£ = a° • Ai,ia° + a° ■ A2,2a° + (a° ■ Ai.sfls + ^2 ■ A2,ia°)P - AqI. (17) 

To get rid of the parity map in the Liouvillean ( !T7|) we define modified hermitian 
Majorana maps 

Sij := ifjalj, a2,j ■= ir/agj-P, a= (0:1,02), (18) 

where r) = 2"i]^."^ai • is a hermitian map obeying the following ( ant i) commutation 
relations 

{al^,fj} = 0, ,[a°^,,r)] = [fi,V] = 0. (19) 

We note that new Majorana maps again satisfy CAR, {a^j,a^^k} = ^pL,v^j,k- The 
Liouvillean takes now a simple quadratic form 

C = a- Aa~ Ao% (20) 

where a structure matrix A is a block matrix 

t: 1:;: ) ■ <^^' 

Incorporating the odd subspace to the Liouvillean enables us to treat the time 
evolution of a general density operator with the formalism of 'third quantization'. Note 
that due to non-normality of the Liouvillean as an operator we have to distinguish 
between left and right vacua 

(1|£ = and £|NESS)=0. (22) 

The left vacuum is simply the identity operator, whereas the right vacuum represents a 
density operator of NESS. 

2.2. Derivation of the Lyapunov equation 

Taking any time dependent solution of quantum Liouville equation p{t) it is convenient 
to study its 2-point correlation function encoded in the 2n x 2n correlation or covariance 
matrix 

Cj,k(t) = ^T^ {wjWkp{t)) - Sj^k, (23) 



Explicit solution of nearly isotropic boundary driven open XY spin 1/2 chain 7 

which may be written compactly in terms of the Majorana maps (llSp as 

C(t) = 2(l|ai®ai|p(t))-l2„. (24) 

Note that in the steady state, since NESS is a generahzed Gaussian sate, all the 
information about NESS can be extracted from C by means of the Wick theorem. 
Expressing the stationarity of the left vacuum ((1| exp(£t) = (1|) and using the trivial 
identitie4§| 

eV^* = i, \p{t))=e%m (25) 

the correlation matrix can be further simplified 

C„fc(t) +5,-fc = 2(l|ai,,ai,fee^*|p(0)) = 2(l|e-^*aije^V^*ai,fce^*|p(0)) (26) 
= 2(l|aij(t)afc(t)|p(0)), 

where a super- Heisenberg picture is defined aj(t) := e~'^*aje'^*. Using quadratic form of 
the Liouvillean fl20l) . we get the Heisenberg equation of motion for the Majorana maps 

,. . X 2 2n 

^^ = [a.,(t),£] = 2 5^^A(.,),(,,)a,,(t), (27) 

M=l 1=1 

where a multi-index (z/, j) is used with z/ = 1, 2 and j = 1, 2, . . . 2n. Differentiating (!26|l 
with respect to time gives 



4 E E( ^(i,,),W)(lKK^)«i,fc(^)|p(0)) (28) 



dt 

At={l,2} 1=1 

+ ^(i,fc),(M,o(l|«ijW<'Wlp(0)) 
Taking into account the identities (assuming / ^ k) 

2{l\au{t)ai4tM0)) = Q^t), 2(l|a2/t)a2,fc(t)|p(0)) = -Q^t), (29) 
2(l|ai,Ki)a2,fe(t)|p(0)) = iQ^t), 2(l|a2/t)ai,fc(t)|p(0)) = iQ,fe(t), 
2(l|a2,fe(t)ai,fe(t)|p(0)) = -i, 2(l|ai,fc(t)a2,fc(t)|p(0)) = i, 

and the definition of A we derive the differential equation for the time-dependent 
correlation matrix 

dC 

-— = -4i[H, C] - 4(M,C + CMj) - 4i(Mi - M;^), (30) 

where M, := Re(M) = (M + M*)/2, Mj := Im(M) = (M - M*)/(2i). Note again 
that the stationary solution ^ = fully determines NESS. The equation for the 
stationary correlation matrix, obtained from (l30!l . can be written compactly using the 
anti-sjTiimetry of Hamiltonian matrix H 

CX^ + XC = Y, (31) 

where new matrices X := 4(iH + Mr) and Y := 4i(Mj'^ — Mi) are defined. This matrix 
equation is known as the continuous Lyapunov equation in control theory and stability 

§ We assume that the Hamiltonian and the couphng operators X^ are time independent. 
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analysis and has a unique solution if and only if the eigenvalues fij of the matrix X 
satisfy [20] 



/3, + /3fc7^0, j,k = l,2,...n. (32) 

We note that although using a general Redfield form of the dissipator the results 
are valid also for the case of Lindblad equation ( 133|) . as the most general Lindblad 
dissipator can be obtained from ([2]) by taking the delta-like bath correlation function 
r^,:y(^) = lfi,uS(t + 0) [H [To]. In the next section we shall apply these theoretical 
considerations precisely in the case of Lindblad equation. 

The Lyapunov equation ( 13T|1 has theoretical and practical implications. On the 
theoretical side, it helps to understand the conditions under which non-equilibrium 
stationary state is unique, it is related to the full spectrum of the structure matrix 
A, and simplifies the analytical calculations, while for practical purposes significantly 
simplifies the numerical implementation of the method of 'third quantization', reduces 
the matrix dimension needed in the computations, and improves the stability of the 
method. 

3. Nearly isotropic XY model 

In the remainder of the paper we consider the Lindblad equation 

^ = -i[H, p] + Y.{2L,pLl - {Lt L„ p}) (33) 

for an open XY spin 1/2 chain with the Hamiltonian 

n-l n 

m=l m=l 

with four Lindblad operators coupling minimally to the ends of the chain 



^l,2 = Vri,2^n ^3,4 = yTJ^S^- (34) 

Using the Jordan- Wigner transformation the XY model is mapped to a fermionic model 
with the Hamiltonian 

n-l ^ -. n 

H = -l^^ i W2mW2m+l H —W2rn-lW2m+2] " i ^ hW2m-lW2m 

m=l m=l 

=: w ■ Hw. (35) 

The equivalent Lindblad operator^ are linear in Majorana fermions 

^1,2 = 2^/^1,2(^1 ±iw^2) =ii,2-w!, (36) 

1 



-^3,4 = 2 V ^l!2(^2n-l ± iW2n) = [3,4 " 'Mi- 
ll The non-local Jordan- Wigner transformation brings a non-local Casimir operator and a phase factor 
in front of £3^4 which however do not alter the Lindblad equation (|33|) . 
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9 



It is convenient to define new coupling constants as T^^ =• ^'^12 ; introducing a small 
parameter e. We also define 



^L,R 



pL,R ^ -p 



L,R 
2 5 



L,R 



L,Ix I L,R, 

rLi ZtZ rto 



(37) 



3.1. Asymptotic analytic solution 

The analytic expression for the correlation function of weakly coupled nearly isotropic 
XY spin 1/2 model can be derived in two steps. Firstly, we make an ansatz for the 
correlation matrix 






ght 



^ 



right 



(38) 



where ^"^ denotes the j-th right eigenvector of X with the corresponding eigenvalue 
/3j, X_^"^ = /3j^j^ ■ Then, by plugging the ansatz into the equation (EI]) we get the 
expression for the coefficients Aj^k 

Aj^k = ^ ^ ^ ^f* ■ Y^f **. (39) 



,T, left * 



Left and right eigenvectors are normalized as ^^ ** ■ ^"^ = Sj^k- 

Secondly, we find approximate eigenvectors of X and determine the coefficients 
Aj fc for nearly isotropic XY spin 1/2 chain with Lindbald reservoirs. Explicit solution is 
found by means of perturbation theory using two main assumptions: (i) small anisotropy 
(7 <C 1), and (ii) weak coupling (6^7). Hence, the matrix X can be decomposed into 
three parts 

X = Xo + 7X^ + eX,. (40) 

The first one is essentially the (isotropic) XX Hamiltonian, the second one is its 
anisotropic part, and the last one comes from the coupling to the environment (4Mi. = 
eX,) 

Xn = 2iha^ !„ - ia^ ® Jr,, 



(41) 



X. 



X. 



/ K^h 










\ 



V . . . 412 / 

where I„ is a ra x n identity matrix and J„ and J'^ are n x n matrices 
/0100...\ /O 1 00 



10 10 
10 1 
10 






-1 





1 





-1 


1 








-1 



(42) 



^. = '-^®tl^ /32. = ±2i(/^-cos(3^)), (43) 



*' n + 1 
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Note that the matrices Xo,X^ are anti-hermitian and real (real and antisymmetric), 
while the matrix X^ is hermitian (real and symmetric). Exact eigenvectors and 
eigenvalues for such system are rather difficult to calculate, therefore, we use the 
perturbation theory to incorporate the effects of the environment and the anisotropy. 
The unperturbed term (Xq) is exactly diagonalizable by the Fourier transformation, 
namely Xq^^- = /S^.i^j with 

where a multi-index (a;,^), a; = 1, 2, j = 1, 2, . . . ,r2 is introduced in order to label 
eigenvalues/eigenvectors, and the plus (+), or minus (— ) sign, in ( HSj) corresponds to 
a; = 1, or a; = 2, respectively, and 

We can neglect the effect of environment to the eigenvectors in the case of very 
small coupling e ^ 7 ^ 1 

%., = ^j + l^,j + 0{e) + 0(f ), (45) 

but we must include the ffist correction to the eigenvalues 

/3-,. = /3°,,+<, + (^(e') + C?(7), (46) 

where, on the contrary, the effects of the anisotropy can be neglected to leading order. 
Reasons for this asymmetry will become clear later on when we will calculate the 
coefficients Aj^. Thus, in case of the eigenvalues we need to calculate only the real 
(dissipative) corrections to ffist order in e, although the imaginary corrections from the 
anti-hermitian part can be much larger (as we have assumed e ^ 7) 

Note that the correction is positive for both, positive and negative imaginary 
unperturbed eigenvalues (energies) Im(/3°j). In contrast to eigenvalues, as already 
noted, important corrections to eigenvectors come from the anti-hermitian part X^ 
(because of the assumption e <^ 7) and are calculated by the usual first order 
perturbation theory 

•,=1 P2,j Pl,j' -,=1 PlJ P2,j' 



with 



i (1 - (-1)^+^') sin (^) sin 



n+1 



Kj,> = ^°*, ■ X^^? ,, = '^\ " '\ '- , "^'\ y^' , . (49) 



sm „r , :{ sm , „. , ,. 

2(n+l) J y 2{n+l) 



Long but straightforward calculation shows that all the matrix elements of X^ 
connecting the vectors with the same polarization (i.e. vectors ^ ■ and ^j/ for 
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j', j = 1, 2 ... n, and a; = 1, 2) vanish. Neglecting environment correction of eigenvectors 
has a convenient consequence, namely, the approximate left and right eigenvectors of 
the matrix X are identical. Note, importantly, that in this approach the anisotropy 
7 must be smaller than 5/3, I7I < 5/3, where 5/3 is the smallest spacing between 
adjacent unperturbed eigenvalues /3° ^ (H3l) which can be estimated as 5/3 ~ ra~^. Thus 
a conservative estimate for the validity of our calculations requires I7I <^ n~^, otherwise 
the degenerate perturbation theory must be used because the difference of the nearest 
eigenvalues can become smaller than the corrections to the eigenvalues. 

The correlation function to the first order in 7 and zeroth order in e may be written 
as a double sum 

^ = Y1Y. A(^,.),K,/)(^, + 7li,,) ® {^',' + 7li',') + 0{e) + 0(f ), (50) 



A(<^j),(a;'j') - e ^ — ^^^ — — ^ — — ^ , [bl) 



with expansion coefficients calculated from equation ( 13911 

(^° J + 7lij) ■ ^'{^',f + 7li',/ 

where Y = eV and 

/ 2K^_ay . . . \ 

... 

Y' 



(52) 



V . . . 2K^a^ ) 

Here an explanation why we used only the environment (hermitian) corrections to 
the eigenvalues is in order. Let us first consider the unperturbed matrix Xq and 
the anisotropy (7) correction X^. Both are real anti-symmetric matrices, hence, the 
eigenvalues of the matrix Xq + 7X^ come in pairs, where one eigenvalue is positive 
imaginary and the other one is equally imaginary negative. Therefore, by considering 
only 7-correction we find singularities in the expression (l50ll . and the condition for 
the uniqueness of NESS (|32|) is not satisfied as well. If we include the symmetric 
(Hermitian), environment part eX^, anti-symmetry of X is broken and the relation ( l50ll 
becomes regular for all j, j' and w, u' . We conclude that the environment (e) corrections 
to eigenvalues are necessary to get a unique NESS. 

Now we will show that indeed only environment corrections to the unperturbed 
eigenvalues are needed in order to get the correct correlation function when the 
inequalities e ^ 7 <^ 1 are satisfied. In this regime the coefficients ( ISTl) are negligible 
unless /3° ,^. = —f3^,j,, which is true for u j^ u' and j = j'. In this case also all 7 
corrections to f3^j come in pairs, as discussed before, and thus add up to zero in the 
denominator of flSTl) . but we have /3K- = P^j^ which remains the only term in the 
denominator. Therefore the small parameter e cancels out of the equation. Hence, 
only first order environment (e) correction to the eigenvalues is needed to calculate 
the "diagonal" coefficients ^(i,j),(2,j), which are the only important contributions to 
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the correlation function, as the rest can be neglected. After simplifying, the largest 
coefficients are found to be constant (i.e. independent of j) 

A(i,,),(2,,) = -A(2,,),(i,,) =: Ao = ^±^ + 0{e) + 0(f ), (53) 

and all the other coefficients vanish in the leading order, i.e. are 0{e^). The correlation 
matrix can now be expressed with a single sum 

n 

C ^ Ao ^ (£^,, ® ^?,, - ^?,, ® ^2,, (54) 

+ l{^l,j ® ^l,j + ^2,j ® ^l,j - ^l,j ® ^2,3 - ^l,j ® ^2,3) 

The error of the previous formula is 0{t) + 0(7^), thus, factors which involve a direct 
product of first order corrections are of higher order and are left out. The dyadic 
products involving zeroth order can be further simplified 

^2 , ® ^? •, - ^? ,• ® ^2 , = ^"^ ® k • ® V" J (55) 

U iJ U iJ 1^ J J J 

as well as the product of the zeroth and first order corrections to the eigenvector^ 

^2j ® ^?,, + ^2,, ® ^l, - ^?,, ® ^L- - ^1,.- ® ^L- (56) 

l—f —J —3 —f 



/-^ /5P _|_ flO 2_7' J_7 2_7 2_7' 



3' 

Summing up the identities fISS]) . f lSB]) and the formula (1511) we get the correlation matrix 
to first order in 7 approximation 

C = Aoa^ ®ln + W ® C" + 0(e) + 0(7'), (57) 

where 

^''■•-^''^''Ew-fw-tr'str (58) 

The final explicit form of the correlation matrix ( 1571) with ( H3|49|l58p is rather 
complicated, but it can be simplified for different regimes. 

First, we observe that analytic properties of the expression Kj^i /{(3^j + Pi^i) in 
(|58|) change at h = h^ := 1. Namely, for \h\ > h^, the expression (l58l) for C^ is analytic 
for any n and can be evaluated asymptotically for large n giving exponential decay of 
the correlator 

CJj, ~ sign(j - f) exp(-|j - j'|/0, n ^ 00, (59) 

where sign(x > 0) = 1, sign(x < 0) = —1, sign(O) = 0, with the correlation length 



^ = l/arcosh(/i) = - 1/ ln(/i - Vh^ - 1) . (60) 

In fig. [H where we show analytically and numerically calculated correlation lengths ^, 
one can clearly see vary good agreement of the explicit formula (|60|) with the numerical 



^ Note a distinction between dyadic product of vectors and Kronecker (tensor) product of matrices, 
both designated with a symbold 0, which should be clear from the type of factors. 
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calculations, and the square-root singularity of the correlation length at the critical field 

h = h„^^ |/l-l|-l/2. 

In the other regime, \h\ < he, the formula ( 157|) can be further simplified in the poles 
of the correlation function (or poles of the coefficients Kjji/{(3lj + /Sij/)), where the 
above non-degenerate perturbation theory fails. This shall be studied in the following 
subsection. 



5 
4 

■jj-. 3 
2 
1 





1.0 



1.5 



2.0 
h 



2.5 



3.0 



Figure 1. Comparison of numerically calculated correlation length ^ (points) with 
the analytical prediction ()60p (solid line). In numerical calculations we used a system 



of size n ~ 20, and with parameters 7 = 10 



= 10- 



Although the size of 



the chain is not yet very large, agreement of numerical coefficients with analytical 
prediction is excellent, except perhaps at the left most point where the finite size 
effect becomes noticable. In this and all other figures we take bath parameters 



d, K2 



1, K 



,R 



d, K2 



The main assumptions that have been made through the calculation are: (i) 
coupling to the environment is the smallest parameter, so the perturbation correction to 
the eigenvectors is negligible, while we must consider the corrections to the eigenvalues, 
(ii) the anisotropy 7 must be smaller than 5/3 so non-degenerate perturbation theory 
around the isotropic case can be applied. 



3.2. Resonances of the correlation function 

In the poles of correlation function the degeneracy of each of the oppositely signed 
unperturbed eigenvalues f3^j, u = 1,2 is twofold. Therefore, we should use the 
degenerate perturbation theory to determine the correlation matrix. As we shall see, the 
sole eigenvectors corresponding to the degenerate eigenvalues determine the dominant 
behavior of the correlation matrix, thus the leading order expression is very simple. The 



)erturbed eigenvalues /3° 


k,l3^i are degenerate (i.e. /3°^ 


— 1^2,1^ f^2,k 


h = hk,i := - 


r f Tck \ / tt/ M 




_ \n + lj \n + lj_ 




some /c, I. 







PI) if 



(61) 
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Let us first assume /c — / is odd, so Kij. 7^ 0. We then diagonalize 'yKij.a^ in the 
degenerate 2x2 subspace to find the proper eigenvectors and eigenvalues: 

C. = ^m^-iy, ~Plk = (3h-lK,k, (62) 

Note again that ^'j' f. and ^ i have orthogonal polarizations fH51) . Let us assume (3i^k > 
{k > I). There is also another degenerate pair with /3i^; = —f3i^k, Ki^k = —Kk,i and 

&,k = ^m?,. - ^Ik), Plk = Ph - lKk,i = -Plk + lK,,k, (63) 

i°, = i=(l?,, + ^°,,), ~Pl = Pi + ^Kk,i = -/??,, - 7^/,fc. 

All the other eigenfunctions and eigenvalues remain the same ^ • = ^ • and /3° ■ = 
/3°j, for i ^ {/c,/}. As in the non-degenerate case, we only need to calculate the real 
(dissipative) corrections to the eigenvalues 

/Si,, = /3ij; forj^jA;,/}, a; = 1,2, (64) 

The coefficients A(i ,,•) (2.j) on the diagonal are calculated from equation (!39|) with new 
eigenvectors and eigenvalues and are almost identical to the non-degenerate case 

A(i,,),(2,,) = Ao = ^^^^, J^{k^} (65) 



and 



where 



sin^ (^) - sin^ (^) 
A(l,fc),(2,fe) = A(l,0,(2,/) = ,-2/^N, -2/^x ^0 = Ao - M (66) 



2sin2(^) 
M := ,, ,, ^"+V/ / ^ Ao■ (67) 

sin^^) + sin^ (;^) 

Diagonal coefficients A(i fc) (2,^) are not constant, therefore, the main contribution to the 
correlation comes from the leading order in the degenerate subspace 

~o ~o ~o ~o ~o ~o ~o ~o ,^„, 

^l,k ® ^2,k - ^2,k ® ^l,k + ^1,1 ® ^2,1 - ^2,1 ® ^1,1 (68) 

= (^^ ®i±l®±l-±^,®±k)■ 
T\le correlation matrix (l50l) has now a simple explicit form 

C,es := C\h=h,, = a^ ® (AoI„ - CL) + 0(7) + 0{e), (69) 

where 

Note that this correction is independent of 7, hence, the correlation in resonance h = h^^i 
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Figure 2. Typical correlation matrix element, C'^/4 3„/4, versus the field strength 
h. In the resonances, h = hk^i, the analytical (orange points) and numerical (black 
circles) correlation element components C^,^ „^,^ are plotted, as the other component 



C 



n/4,3n/4 



is orders of magnitude smaller. For other values of the field we compare 
analytical (orange dashed line) and numerical (black solid line) correlation element 
components C^,^ 3n/4' whereas here the component C^,^ „^,^ is considerably smaller 



(see fig.[3|). In both cases the agreement is excellent. System parameters are: 
7 = 10-4^ e = 10~6. 



20, 



( 16T]) scales as 1/n - this is hidden in the normahzation of ip . ( H3|) . It is also interesting 
to note that in the leading-order at the resonance the correlation matrix ( !70|l is a simple 
combination of two 'sine-waves' and is thus, for large n, very much reminiscent of 
an eigen-mode of a square drum. This is why we choose to call this phenomenon a 
correlation resonance. 

In order to obtain the next-to-leading first order 0(7) correction in the resonances 
we successively apply the degenerate and non-degenerate perturbation theory. The 
corrected eigenvectors are then manipulated as in the previous cases resulting in a 
rather complicated formula 



where 






2iAo Yl 

^<j,j'<n 






—J 



(71) 



(72) 



+ i{Ao-6A) J2 
- i(Ao - SA) Y. 



K^ 



k,j 



Ph + Ph^-^^-''-'"'-^ 



'l,j 






Ph + Ph -^ 



{±. 0tp,-^,0^ 



J-l Z.I 



—3 
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where the coefficient G is obtained from the restriction of X in the degenerate subspace. 
In the case when Ki^k = 0, which is true for even-even or odd-odd k, I, the degenerate 
eigenvectors are not directly coupled, so the general result is only slightly changed 
with respect to the non-degenerate case ( !57f58l) . namely one only needs to exclude the 
degenerate pair from the double summation 

Cres = Ao^" ® In + 7^" ® C^,, + 0{e) + 0(7^), (73) 

where 

CL = 2iAo E .0 fao t,®tf (74) 

In fig. |2] we compare the analytical predictions fl57|) off the resonances and (|69|) on 
the resonances with numerical simulations and find excellent agreement in the entire 
long range order region \h\ < 1. In fig. [3] we zoom in a narrow window around a single 
resonance and see very clearly that the resonance width scales as 

5h K, 7, (75) 

i.e. the formula fl57|) is essentially applicable for \h — /ifc^l > 5h, whereas the formulas 
(EHl EH EHD hold ioi\h- hk,i\ < 5h. 

Before closing, several interesting remarks are in order. Firstly, we note that in 
resonances the correlation function C^ ( !70|) does not depend on the anisotropy 7, 
however the use of perturbative arguments indicates that this is true only for small 
7, i.e. I7I < Sp. As noted before, for larger but still small anisotropics (1/n^ ^7^1) 
degenerate perturbation theory must be used as some pair of unperturbed eigenvalues 
/3° ■ may become nearly degenerate. In other words, for fixed small anisotropy, our 
results are quantitatively valid only up to a maximal size nmax ~ ItI^^^^- However, 
numerical calculations (see fig. Hj) indicate that this estimate is too conservative, in fact 
7 can be as large as ~ 1/n for our results to remain valued. This means that, in practice, 
6/3 can be estimated as an average spacing between eigenvalues ( l43l) . rather than the 
more conservative minimal spacing. Nevertheless, the fact that we cannot extend our 
calculation to thermodynamic limit n — )■ 00 for any fixed system parameters implies 
that our method can not be used to describe the critical closing of Liouvilean spectral 
gap behavior as observed in [9l [TO] . 

Secondly, we stress again an interesting fact that for small reservoir coupling 
strength e our leading order perturbative results do not depend on e, meaning for 
example that the value of the spin-spin correlation at long distance in the regime \h\ < he 
is insensitive to e even for infinitesimal coupling, but there is, however, dependence 
on the properties of the baths (e.g. effective temperature, chemical potential, etc) as 
encoded in parameter Aq. We note that in the models of symmetric out-of-equilibrium 
driving at infinite temperature (e.g. [121 1131 [21]), where formally Aq = 0, one obtains 
a different result: C = 0{e). With increasing e, our results remain valid - due to our 
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perturbative assumptions - up to e* ~ 7. This is illustrated with simple numerical 
calculations in fig. [51 

The last point we make concerns the total number N{n) of resonant values of the 
magnetic field hk^i for a given chain length n. A simple inspection of the formula (16T|) 
with the constraint that k — I should be odd, results in an expression for the resonance 
counting number 

N,,,{2n) = N,,,{2n-l)=n{n-l). (76) 

If we consider only the values of the field with a fixed, say positive sign, h > 0, then the 
number of resonances is N{n)/2. This amounts to, for example, in total of 45 resonances 
shown in fig. [2] for n = 20. 
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Figure 3. Detailed plot of the resonance around /ii,2 = (cos^(^^) + cos (^^)/2 for 
a system with n — 20, e = 10~^ and two values of the anisotropy, 7 — 0.5 x lO^'' (black 
curves) and twice larger 7 = 10^^ (orange curves). The dashed curves show numerically 
obtained C^u 3^/4 • The solid curves show C,^/4 3„ /4, where the thick curves correspond 
to numerical results and the thin ones to analytical calculation, and the blue, and gray, 
circles corresponds to analytically predicted C^^ 3„/47 ^^d C^/4 3„/4, at the resonance 
field, respectively. The resonance width grows linear with the anisotropy 7, also the 
agreement between the (non-resonant) analytically and numerically obtained values of 
a^ projection (C^) is good for \h — hi^2\ > 7- 



4. Discussion and conclusion 



The present article contains two main results. Firstly, using the quantization in the 
Fock space of operators ('third quantization') we have derived manifestly bilinear form 
of the many-body Liouvillean and expressed the equation of motion for the 2-point 
correlation function in terms of the continuous Lyapunov equation. Secondly, we have 
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Figure 4. Comparison between analytical expression (thin lines) and numerical 
calculations (thick curves) for the correlators C^u 3„/4 (sohd lines) and C^u 3„/4 
(dotted lines) versus the anisotropy parameter 7. Results for n = 20 (black), n = 40 
(blue), n = 80 (orange) are plotted. Note that results show clearly that perturbative 
analytics accurately describes the numerics up to 7 « 1/n. Horizontal thin lines are 
calculated from the equation ([70|) and the diagonal (dashed) thin line from the formula 
(El. 
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Figure 5. Comparison between analytical expression (thin horizontal lines) and 
numerical calculations (thick curves) for the correlator C^m 2,nli "^^rsus the reservoir 
coupling strength e, and for two different values of anisotropy 7 = 10~^, 10~^ (black, 
orange, respectively) indicated with vertical dashed lines. Solid lines refer to C^ 
component and dashed lines to C^ component of the correlation matrx. 
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presented explicit analytical result for the steady state correlations in boundary driven 
open XY spin 1/2 chain in the Lindblad formulation. The two regimes of short range 
and long range correlations have been clearly identified, with the non-equilibrium phase 
transition occurring at the critical value of transverse magnetic field /ic = 1- In the short 
range regime, exponent of correlation decay has been obtained analytically, whereas in 
the regime of long range order, a particular resonant spikes of the correlation response 
of the system have been identified, where the values of the correlator (scaling with the 
inverse of the chain length n) have been calculated in a simple closed form. 

We note that our analysis could perhaps be extended to the case of large (or better 
to say, non-small) anisotropies, where the continuous limit (for n — ?> oo) of the Lyapunov 
equation of the XY model becomes a non-homogeneous Helmholtz equation driven by 
two point-like sources representing the two Lindblad reservoirs. The interpretation of 
the non-equilibrium phase transition in this picture becomes rather obvious, namely 
long-range order corresponds to positive energies and real wave numbers, whereas short 
range phase correspond to negative energies and imaginary (evanescent) wavenumbers. 
The correlation resonances can then be exactly identified with the eigenstates of a 
square-shaped quantum billiard (with an appropriate symmetry/boundary contition). 
Nevertheless, the exact analytical results for the correspondence between the open 
(strongly anisotropic) XY chain and the driven Helmholtz equation remains a work 
in progress. 

We hope that our results, which are nontrivial though surprisingly simple, may be 
relevant and observable in real laboratory experiments. 
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